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We study the rapidity window dependences of higher order cumulants of conserved charges ob¬ 
served in relativistic heavy ion collisions. The time evolution and the rapidity window dependence of 
the non-Gaussian fluctuations are described by the diffusion master equation. Analytic formulas for 
the time evolution of cumulants in a rapidity window are obtained for arbitrary initial conditions. 

We discuss that the rapidity window dependences of the non-Gaussian cumulants have characteristic 
structures reflecting the non-equilibrium property of fluctuations, which can be observed in rela¬ 
tivistic heavy ion collisions with the present detectors. It is argued that various information on the 
thermal and transport properties of the hot medium can be revealed experimentally by the study 
of the rapidity window dependences, especially by the combined use, of the higher order cumulants. 
Formulas of higher order cumulants for a probability distribution composed of sub-probabilities, 
which are useful for various studies of non-Gaussian cumulants, are also presented. 
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I. INTRODUCTION 


Bulk fluctuations of conserved charges are believed to be useful observables for the study of thermal property of the 
primordial medium created in relativistic heavy ion collisions. Active experimental analyses of fluctuation observables 
have been performed in the beam-energy scan (BES) program [I| at the Relativistic Heavy Ion Collider (RHIC) M, 
as well as the Large Hadron Collider (LHC) These experiments have been observed not only the second order but 
also higher order cumulants characterizing non-Gaussianity of the fluctuations by the event-by-event analysis. Various 
theoretical suggestions have been made to utilize fluctuation observables as signals of the deconfinement transition 
and the QCD critical point din. The cumulants of conserved charges in the equilibrated medium are also observable 
in lattice QCD Monte Carlo simulations, and the numerical study of various cumulants have been carried out actively 


One of the most notable features of the measurement of fluctuations in heavy ion collisions is that higher order 
cumulants of particle numbers in a rapidity window have been observed with moderate statistics up to fourth order. 
Such measurements are possible because the systems observed in these experiments are not large; the particle number 
measured in one event is at most of order 10 ^ although the local equilibration is expected to realize for some 

period in the time evolution. Experimental detectors are designed to observe and identify almost all charged particles 
entering the detector. They thus can provide full counting statistics (sTj l of the event-by-event distribution of the 
particle numbers up to the error arising from the acceptance, efficiency and particle miss identifications [T^ . l32l - [3^ 
of the detectors. 

Because the experimental detectors in heavy ion collisions count particle numbers which arrive at the detectors, 
the event-by-event analysis can observe fluctuations in the final state of the collision events. If the fluctuations are 
fully equilibrated in the final state, they can only tell us information on the thermodynamics at this time, which 
is in the hadronic medium and well described by the hadron resonance gas (HRG) model [T^ |. The experimental 
results, however, suggest that the observed fluctuations are not the one of the equilibrated hadronic medium. Eirst, 
the experimental results on the cumulants and their ratios have statistically-significant deviation from the values 
in the HRG model ii- In the HRG model, the fluctuations of net particle numbers, including net-baryon and 
net-electric charge which are conserved charges, are given by the Skellam distribution to a good approximation [l3l| . 
The cumulants observed at RHIG and LHC show deviations from these Skellam values. The second observation is 
concerned with the rapidity window dependences of the cumulants. In a thermal medium in grand canonical ensemble, 
cumulants of conserved charges are extensive variables and proportional to volume. In heavy ion collisions, fluctuations 
of particle number in a given pseudo-rapidity window, Arj, are observed. Assuming the Bjorken expansion and that 
the pseudo-rapidity window can be used as a proxy of the one of the coordinate-space rapidity. Ay [s^, the spatial 
volume to count the particle number is proportional to Arj. The cumulants, therefore, should be proportional to A 77 if 
the fluctuations are fully equilibrated. The net-electric charge fluctuation {{SNq)'^) observed by ALICE collaboration 
0, however, clearly shows that this proportionality is violated. A similar violation is also indicated in the higher order 
cumulants of net-proton number observed by STAR collaboration [^ . It is this non-equilibrium nature of fluctuations 
which enables us to explore the primordial thermodynamics using fluctuation observables [^, . 


For conserved charges, the off-equilibrium property of fluctuations in the final state is not surprising, because the 
modification of local density of a conserved charge is proceeded only by the diffusive processes, which is typically slow; 
because of the conservation law, the time scale of the diffusive process can become arbitrary slow as the spatial scale 
becomes larger. In fact, earlier studies suggested that the time scale of the diffusive process for conserved charges is so 
slow that the thermal fluctuations generated in the deconfined medium survive until the final state when Arj is taken 
sufficiently large [1, |^, Is^, . Because the effect of the diffusion with a fixed Ay is determined by the transport of the 
charges, the conservation law also suggests that one can investigate the transport property of the hot medium, such 
as the magnitude of the diffusion constant, experimentally using the Ay dependence of conserved-charge fluctuations 


In Refs. Hi 113,131, the time evolution of fluctuations are considered only for Gaussian fluctuations. Because 
the heavy ion collisions can measure the higher order cumulants, the Ay dependences of not only Gaussian but also 
non-Gaussian fluctuations should encode more information on the thermal and transport properties of the medium 
M- To extract such information, however, an appropriate description of the non-equilibrium time evolution of the 
non-Gaussian fluctuations is required. 

When one investigates the diffusive process of non-Gaussian fluctuations in heavy ion collisions, it is important to 
keep in mind that the cumulants observed so far are not far from the ones of the Skellam distribution ill , which are 
the equilibrated values in the HRG model. A possible interpretation of these results is that the signals in cumulants 
are developed in the primordial stage, but they are blurred by the diffusive process in the hadronic stage. In fact, the 
experimental result in Ref. [3 in dicates that the fluctuations increase toward the equilibrated value in the hadronic 
medium due to the diffusion [IqI [13, HI ■ If this is the case, it is important to describe the approach of the cumulants 
toward the hadronic value in the analysis of the diffusive process in the hadronic medium. For this purpose, the model 
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FIG. 1. System described by the diffusion master equation (DME) Eq. ([Tjl. 


for the diffusive process should have nonzero higher order cumulants in equilibrium which are consistent with those 
in the hadronic medium. 


In Ref. to describe the time evolution of non-Gaussian fluctuations the diffusion master equation (DME) is 
employed. As discussed in Ref. , in this model the fluctuations in equilibrium are given by the Poisson or Skellam 
distribution as a consequence of the discrete nature of the particle number. The time evolutions of the second and 
fourth order cumulants of net particle numbers are studied in this model. It is found that the fourth-order cumulant 
can show characteristic behaviors as a function of Ary and tends to be suppressed compared with the second order one, 
if the thermal fluctuations in the primordial medium are suppressed. These results can be confirmed experimentally 
with the present detectors. The analysis is later extended to finite volume systems 2l| to investigate the effect of the 
global charge conservation. 


The purpose of the present paper is to elaborate on detailed discussions skipped in Ref. [T^. In the present study, 
we show the analytic procedure for obtaining the time evolution of cumulants in the DME with arbitrary initial 
conditions. We also extend the discussion on the Ary dependence of cumulants in Ref. [l^ to more general cases. The 
Ary dependence of the third order cumulant is considered in addition to second and fourth order ones discussed in 
Ref. [1^. The initial condition is also extended to more general cases. Using these results, we pursue the possibility 
to extract various information on the thermal and transport properties of the hot medium from the non-Gaussian 
cumulants of net-baryon number and net-electric charge in heavy ion collisions. 


This paper is organized as follows. In Sec. |nl we present the analytic solution for the time evolution of cumulants 
in the DME. The result is then extended in Sec. m to the multi-particle systems in order to describe the cumulants 
of net particle numbers. In Sec. m we apply these results to the description of Arj dependences of non-Gaussian 
cumulants in heavy ion collisions. Some discussions and a short summary is given in Sec. [V] In Appendix we 
derive formulas of cumulants for a probability given by a superposition of sub-probabilities, which are used in Sec. [TTl 
to obtain the solution of the DME for arbitrary initial conditions. In Appendix |B1 we consider a simple stochastic 
model for a chemical reaction. This appendix will help readers who are not familiar with stochastic models like the 
DME to understand Sec. HIl 


II. DIFFUSION MASTER EQUATION 


To describe the diffusion of non-Gaussian fluctuations of the conserved charges in heavy ion collisions, we use the 
diffusion master equation (DME) [I^. After introducing the DME and describing its general property in Sec. Ill Al 
we derive the analytic form of the time evolution of the cumulants for fixed initial conditions in Secs. Ill 51 - III FI The 
solution is then extended to general initial conditions in Sec. IIIGI 


A. Diffusion master equation 


We consider a single species of classical Brownian particles moving in one dimensional space. To describe microscopic 
states of this system, we divide the spatial coordinate into discrete cells with an equal length a, as in Fig. [TJ We label 
each cell by an integer m and denote the number of particles in each cell by Um- A microscopic state is then specified 
by the vector n = (• • • , rim-i, n-m, n-m-i-i, • • •) representing particle numbers in all cells. One can also introduce the 
probability distribution function P{n,t) that each cell contains rim particles at time t. 

Next, we assume that each particle moves adjacent cells with a probability 7 (t) per unit time. The probability 




















4 


P(n, t) then follows the first-order differential equation 

d 

^ “h Gjyi^i G'lji—l') “h P{tI Gjyi^i p 

m 

-2n„P(n,t)], (1) 

where = (0, • • • , 0,1, 0, • • • , 0) is the vector that all components are zero except for the mth-one which takes unity. 
After solving Eq. o, we take the continuum limit defined by a —>■ 0 to recover the continuity of the spatial coordinate. 
Although the generalization of the model to a multi-dimensional system is straightforward, in this study we limit our 
attention to the one-dimensional problem to avoid unnecessary complexity. Note that the one-dimensional model is 
sufficient for the description of the diffusive process along the rapidity direction in heavy ion collisions. 

As we will see later, after taking the continuum limit the time evolution of the average particle number density 
n(a;, t) (the exact definition will be given in Sec. IIIDj) in the DME obeys the standard diffusion equation 

d ^2 

—n(a:, t) = D{t)—n{x, t), (2) 

with the diffusion constant D{t) = 'y(t)a‘^. It is also shown [l^ that the time evolution of the Gaussian fluctuations 
of n{x, i) in this limit is consistent with the stochastic diffusion equation 

-^n{x,t) = — (D{t)—n{x,t)+£,{x,t)] , (3) 

where (,{x,t) represents the stochastic term. Here, it is assumed that ^{x,t) is local with respect to x and t, i.e. 
(^(xi, ti)^(a: 2 , ^ 2 )) is proportional to i5(xi — X2)5{ti — t 2 ). In this case the property of ^(x, t) is completely determined 
by the fluctuation dissipation relation [40l |. 

The DME, on the other hand, has a property that the higher order cumulants take nonzero values in equilibrium 
defined by t —>■ 00 limit. This property is contrasted with Eq. m in which the fluctuation of n(x, t) becomes Gaussian 
and all the cumulants of n(x,t) higher than the second order vanish in equilibrium for n independent D{t) [T^. The 
higher order cumulants of the particle number in some range of the spatial coordinate in the DME in equilibrium are 
given by the Poissonian ones [l9l | . This property is easily understood from the fact that in equilibrium each Brownian 
particle exists any place with an equal probability and the individual particles are uncorrelated. When the model is 
extended to treat the net-particle number in Sec. EH the distribution is given by the Skellam one. This property is 
suitable to model the time evolution of higher order cumulants in heavy ion collisions. 


B. Factorial generating function for fixed initial condition 


In the following, we solve the DME Eq. O- For readers who are not familiar with stochastic equations like Eq. o, 
we recommend them to read Appendix iBl before this section. 

In the following hve subsections, we first solve Eq. m for the fixed initial condition. 


= (4) 

m 

and present t dependence of cumulants of P{n,t). In Eq. (|4]), the particle number in each cell, Um, are fixed to Mm 
without fluctuations at t = 0. The analysis is later generalized to initial conditions having fluctuations in Sec . Ill Gl 
To solve Eq. © with the initial condition Eq. ®, it is convenient to use the factorial generating function 

(5) 


( 6 ) 


Gf (s, ^)= XI (n 

n ^ m ' 

where the sum is taken over all possible combinations of rim- Substituting Eq. (O in Eq. o, one obtains 

d d 

(s, f) = y{t) ^ -f Sm—i)~T^ Gf (s, t). 


dSr, 


Because Eq. is a first-order partial differential equation, this equation is solved by the method of characteristics. 
In this method, we use the fact that the solution of Eq. m satisfies 


|Gf(s^(t),t) = 0, 


( 7 ) 
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where the characteristic line (t) is the solution of the characteristic equations 

(j 


dt 


= -7(0(Sm+l - 2s^ + Sm-l)- 


Equation (jS]) is easily solved in Fourier space. Assuming the periodic boundary condition and defining 

1 

N ■ 


~ /V X! 


I'KijmjN 


with TV being the total number of cells, the characteristic equations become 

d 

dt ^ 

with 




27T 


where the last nearly-equality is satisfied for j /N <C 1. The solution of Eq. (jTU)) is 


rj (t) = Tj (0) exp[ / dt'ujj {t')]. 
Jo 


When 7 (t) does not have t dependence, Eq. (fT^ has a simple form 

r|(t) = 

The initial condition Eq. o is transferred to the factorial generating function G{{s,t) as 

M„, 

G({s,0) = = n ' 

m m \ j 

Because of Eq. 0, the solution of Eq. ([H) is obtained by substituting Eq. (IT^ in Eq. (HI which yields 


M„ 


=n E 


T-jexp[— J dt'ujj (t')] 


m \ J 

The corresponding factorial cumulant generating function is given by 


K{{s, t) = log G{{s, t) Mm log ^ 


^2mjm/N 


rie"""-'""'*’ exp[— / dt'u)j{t')\. 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 


(15) 


(16) 


C. Cumulants and factorial cumulants 

Using Eq. (HU, the factorial cumulants of Um are given by 

(^mi'^m2 ' ' ' ^m/)fc — 


d^Ki 


m2 ' ' ' ^^mi 


3 = 1 


(17) 


where s — 1 means that = 1 for all m. Note that this condition corresponds to tq = 1 and rj = 0 for j / 0. Using 
this relation, factorial cumulants in the Fourier space are calculated to be 


{jlk^Tlk2 * ‘ ‘ ^ki)ic — ^ ^ ^ 

mi,-" , 771 ; 

= E ^ 

mi,-" , 771 ; 

d^Ki 


-2-!Ti{kimi + ---kimi)/N 1 ~ . . . \, 

X'^TTli '^771;/IC 

d^Ki 


— 27ri(fei77iiH— kimi) /N 


f)c . . . c 


drki ■ ■ ■ drki 


(18) 
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where the Fourier transform of the particle number is defined by 


1 


X k 

Substituting Eq. (HU) in Eq. (US, one obtains the factorial cumulants up to fourth order as 


(fife)fc =Mfc exp[- / dt'ujkit')], 

Jo 

(fifei fife 2 )fc — ^k\+k 2 i.^ki (^ ) “t“ ^k 2 id' ))] 5 

Jo 

{nki n^k2^k3 ) fc — 2 A^fe^_|_fe 2 +fe 3 exp[ / dt i^ki (t ) ^k^ii ))]? 

JQ 

{nkink^nkifiki)fc = - 6Mki+k2+k3+kiexp[- / dt'{iOk^it') + ujk2it') + uJtsit') + ujk^{t'))], 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 
(23) 


with Mfe = Although we do not show the explicit forms of cumulants higher than fourth order, 

they can be obtained straightforwardly with a similar manipulation. 

Cumulants of fife are obtained by using cumulant generating function 


as 


K{6,t) = K{is,t)\^.^^ei, 
. _ d^K 

\^mi^m2 ^mi/c — 


56>i • • • dOi 


0^0 


Using the relations 


dSm _ ^ d Sm _ ^ ^ 

— ^m,m'•) mi 7712 5 

^^7712 


de. 


(24) 


(25) 


(26) 


and so forth, one finds that the cumulants of fife are related to factorial cumulants Eqs. (Pn ]) - (P I) as 


{^k'}c -(^fc)fc7 

{riki^k'^c — ('hfcj^'hfc2)fc H“ (hfcj^-|-fc2)fc; 

'h/i ;2 h/eg)c ^A; 2 ^^ 3 )fc H“ {^^k\^k 2 -\-k-^ic (comb.) 

+ (h/ei + fc2+fe3)fc) (29) 

(fifei fife2fife3fife4)c =(fifcififc2^fe3^fe4)fc + (fifci fifes ^/c3+fe4)fc + (comb.) 

+ (fifei+fe 2 fife 3 +fe 4 )fc + (comb.) + {nkink^+k^+ki)^ + (comb.) 

+ (fifel+fc2+fe3+fe4)fci (30) 

where (comb.) means the sum for all possible combinations for subscripts of fife; for example, 

(fifeififes+fe3)fc “t“ (comb.) — (fifeifife2+fe3)fc (fifesfifes+fei)fc (fifes^fei+fe 2 )fc- (31) 


D. Continuum limit 

Next, we take the continuum limit, a 0. First, we denote the spatial coordinate of the mth cell as cc = ma. In 
the continuum limit, using the particle number density defined by n(x) = rimja, the probability distribution function 
in Eq. (HD is promoted to the functional P[nix),t] of n{x). Note, however, that this notation is conceptual; in actual 
applications, the functional P[nir]),T] is always understood as the limit of the function P{'n,T) with small but finite 
a. The Fourier transform of nix) is n{p) = f dxe~'^J‘^n{x) with p = 2'Kk/Na. One can easily verify that the factorial 
cumulants of n{p) in the continuum notation are obtained by replacing fife —>■ nip) in Eqs. (l20l) - (|2^ with 

ojpit) = 7(<)aV. 


(32) 
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We require that the deterministic part of Eq. m obeys the solution of the diffusion equation Eq. ([2]) in the continuum 
limit. The solution of Eq. m in Fourier space with the initial condition n{x, 0) = M (x) is 

h{p,t) = M{p) exp[— f dt'D{t')p‘^]. (33) 

Jo 

By comparing Eq. (1331) with Eq. (12(111 . one finds that these equations give the same result when one takes 

D{t) =. (34) 

The a —>■ 0 limit therefore has to be taken with Eq. (IMl) . Similarly, one can show that the time evolution of the 
second order cumulant in the DME, Eq. (ED, in this limit agrees with the solution of the stochastic diffusion equation 

Eq. dSD [H. 


E. Particle number in a volume 


Next, we investigate the cumulants of the particle number in a range with length A, 

^A/2 


Q{A,t) = / dxn{x,t). 

J-A/2 


(35) 


To represent the cumulants of Q, we first note that the diffusion distance of each Brownian particle in the model is 
given by 


d{t) = -^2 dt'D{t'), 


(36) 


from Eq. (|33)) [21|, Using this quantity having the dimension of length, one can define a dimensionless variable 

d{t) 


X = 


A 


(37) 


The t and A dependences of (((5(A, t))”)c/A then are described only through this dimensionless variable, because 
(((5(A, t))"’)c/A is dimensionless and Eq. (1571) is the natural dimensionless variable which can be created from the 
combination of A, t and D{t). 

Dependences of the first order cumulant on A and t is calculated to be 

{Q{A,t))c= / dx{n{x))c = / dx :r^e*P^(n(p))c 
J-A/2 J-A/2 J 

/•^/2 r 

= dx y-e*P“M(p)e-“^‘ 

J_A/2 J 27r 
A/2 


= J dzM{z)Ix{z/A). 


^ip(x-z)^-Dtp 


(38) 


Similarly, the second order cumulant of Q{A,t) is calculated to be 

((Q(A,t))^)c = / dxidx 2 [ 

J-A/2 J 27r)^ 


-A/2 

fA/2 


(27r)2 

dxidx2 f 

J-A/2 J (271-)^ 

X M{pi +p2)(e““’’i+*’2‘ - e“(‘^Pi+“P2)*) 

= f dzM{z) f ' dxidx 2 f 

J J-A/2 J (27r)^ 

X _ g-Dtipl+pl)^ 

= J dzMiz) {Ix{z/A) - Ix{z/Af) . 


(39) 






To obtain the last line, it is convenient to use a relation 


+P2) = S{X, - X2) [ 


( 2 ^) 


(40) 


In the last line we have defined 


/v«) = f'\ j = 1 (^ert 




y/2X 


y/2X 


(41) 


where erf (a;) = (2/y^) fjdte * is the error function. Some properties of Eq. (j41l) is summarized in appendix [Cl 
Repeating similar manipulations, one finds that higher order cumulants of Q{A,t) up to fourth order are given by 


((Q(A,t))")e= / dzM{z)H„{z) 


(42) 


with 


H,{z) =Ix{z/A), (43) 

H2{z)=Ix{z/A)-Ix{z/A)^, (44) 

H^iz) =Ixiz/A) - iIx{z/Af + 2Ix{z/A)\ (45) 

Hi{z) =Ix{z/A) - llx^zjAf + 12/x(z/A) 3 - 6Jx(^/A)^ (46) 


F. Uniform initial condition 

As a special case of the above result, it is instructive to consider the solution of the DME for a fixed initial condition 
with a uniform density, n(x, 0) = M(x) = M. With this initial condition, the average density trivially takes a constant 
value for all t. The fluctuations, however, have nontrivial t dependence even for this simple case, because they increase 
toward the equilibrated value from zero as t becomes larger. 

Substituting the initial condition M{x) = M in Eq. (1421) one obtains 


(Q(A,t))e=MA, (47) 

((Q(A,t))2),=MA(l-F2(A)), (48) 

((Q(A, t)f), =MA (1 - 3E2(X) + 2F3(A)) , (49) 

((Q(A, t))4), =MA (1 - 7E2(X) + 12E3(A) - 6F4(A)), (50) 

with 

Fn{X) = j dz[Ix{zlA)\^. (51) 


Some properties of Fn{X) are given in Appendix [C] Note that the A and t dependences of the cumulants are encoded 
through A; as in Eq. (EZl), A = 0 corresponds to t = 0 while A goes to infinity for t —>■ oo with a nonzero A. 

Eq. (lT7l) shows that the average (Q) = {Q)c = MA does not have t dependence. On the other hand, the second 
and higher order cumulants, Eqs. (H51) - (|5(I)) . are t dependent. As discussed in Appendix [Cl Fn{X) for n > 2 are 
monotonically decreasing functions of A with Fn{0) = 1 and lim^^oo E'„(A) = 0. Erom these properties, one finds 
that Eqs. (1481) - (I5II1) vanish at t = 0 as are required by the fixed initial condition, while they approach the Poissonian 
value {{Q{A,t))'^)c = MA in the large t limit. 

In Eig. [21 we show the second, third and fourth order cumulants normalized by their equilibrated value, MA, as 
functions of A and 1/A. Because A is a monotonically increasing function of t with fixed A, the left panel can be seen 
as the time evolution of {{Q{A, t))'^)c. The panel shows that the second and third order cumulants are monotonically 
increasing functions. The fourth order one, on the other hand, has a non monotonic time evolution; it first becomes 
negative, and then increases toward the equilibrated value. It is worth emphasizing that the fourth order cumulant 
can become negative owing to the non-equilibrium effect. One also finds from the figure that the approach of the 
cumulant to the equilibrated value is slower for higher order. An intuitive interpretation of this tendency is discussed 
in Appendix |BJ 
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FIG. 2. Dependences of the cumulants of Q(A,r) on X = d{t)/A (left) and 1/X (right). 


In the right panel of Fig. [51 the same quantities are shown as functions of 1/X = A/d{t). With fixed t, this plot 
can be regarded as the A dependence of {{Q{A,t))^)c/MA. The figure shows that the equilibration of (((5(A, f))”)c 
is realized in the A —>■ 0 limit, while {{Q{A,t))'^)c approaches the initial value as A becomes larger. In other words, 
the approach to the equilibrated value is faster for smaller A. We note that such a behavior of second-order cumulant 
of net-electric charge fluctuation is observed by the ALICE collaboration Q, as will be discussed in more detail in 
Sec. HVCl 


G. General solution 


So far, we have considered the DME Eq. o for the fixed initial conditions, n(x,0) = M(x) without fluctuations. 
Now, we extend these results to general initial conditions having fluctuations. This procedure is nicely carried out by 
a superposition of the previous results and using the formula of cumulants given in appendix [^ 

Let us first denote the probability distribution in the initial condition at f = 0 as P[n(a;),0] = F'[n(a;)]. Because 
Eq. m is a linear equation, probability distribution P[n(x),t] with the initial distribution P[n(a;),0] = F'[n(a;)] is 
then given by 


P[n{x),t]= ^ F[M{x)]PM{x)[n{x),t], 


(52) 


where PM(x)[n'ix)yt] is the probability distribution function at time t with the fixed initial condition with n(z, 0) = 
M (x), and the sum runs over the function space of M (x). 

Eq. (1551) shows that P[n{x),t] is given by the superposition of probabilities Pm (x)['n‘{x), t] with an weight F[M{x)]. 
In addition, the cumulants of PM{x)[n{x),t], given in Eq. (H5]) . are linear with respect to M{x). Therefore, the 
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cumulants of Eq. (l52)) is obtained by substituting Eq. (1421) to Eqs. (IA20p - (IA23I) . which results in 

{Q{A,t))c = J dz{M{z))c,oHi(z), (53) 

{{Q{A,t)y)c = J dzidz2{M{zi)M{z2))c,oHi{zi)Hi{z2) 

+ J dz{M{z))cfiH2{z), (54) 

{{Q{A,t)f)c = J dzidz2dz3{M{zi)M{z2)M{z3))c,oHi{zi)Hi{z2)Hi{z3) 

+ 3 j dzidz2{M{zi)M{z2))cfiHi{zi)H2{z2) 

+ J dz{M{z))c,oH3iz), (55) 

r ^ 

{{Q{A,t))*)c = / dzidz2dz3dz4{M{zi)M{z2)M{z3)M{z4))c,oY[Hi{zi) 

i=l 

+ 6j dzidz2dz3{M{zi)M{z2)M{z3))cfiHi{zi)Hi{z2)H2{z3) 

+ J dzidz2{M{zi)M{z2))c,Q{3H2izi)H2iz2) + AHi{zi)H3{z2)} 

+ [ dz{M{z))c,oH 4 {z). (56) 


where {M{zi)M{z 2 ) ■ ■ ■ )c,o represents the cumulants of the initial distribution F\M{z)]. Note that the above results 
up to second order are consistent with the solution of the stochastic diffusion equation. 


H. Initial condition satisfying uniformity and locality 


Using Eqs. (1531) - (1561) one can obtain cumulants of Q{A,t) with an arbitrary initial condition having fluctuations. 
We next constrain the initial condition so that the result becomes physically apparent. 

In an equilibrated medium and in a macroscopic scale, density fluctuations are local, i.e. the cumulants of the 
particle density is given by 


(n(a;i)n(a;2) • • • n{xi))c =[n']c5(a:i - ^2) ■ ■ • 5{xi - Xi), 
where the coefficients in Eq. dSZD defined by 


(57) 


K], = (Q")./A, (58) 

represent the generalized susceptibility. This is because the correlation length of the density would be at most the 
same order as the microscopic scales such as the mean free path without long range correlations. The correlation 
between macroscopically separated points, therefore, must vanishes. It is notable that Eq. (l57l) is also satisfied in the 
classical free gas in grand canonical ensemble. When, and only when, the cumulants satisfy the condition Eq. (113, 
{Q'^)c are extensive variables and proportional to A, 

I'A 

((Q(A,t))^)c = / dxi - • dxn{n{xi)n{x2) • • • n{xi))c,o = A[n^]c. (59) 

It is reasonable from Eq. (l59l) to refer to [n^\c as the generalized susceptibility. Note that the translational symmetry 
of the system is assumed in Eq. (113; extension of Eq. (1571) to the finite volume system is addressed in Ref. [2l| . 






11 


For initial conditions satisfying Eq. (153, the cumulants Eqs. (ESI) - Ei) are reduced to be 

(Q(A,t))e=A[M]„ (60) 

((Q(A, t)f), =A[M], {1 - F^iX)} + A[M\F2{X), (61) 

((Q(A, t)f), =A[M], {1 - 3F2(A) + 2F3{X)} + :iA[M% {F2{X) - F^iX)} 

+ A[M%F3iX), (62) 

((Q(A,t))4), =A[M], {1 - 7F2(A) + 12F3(A) - 6F4(X)} 

+ A[M% {7F2iX) - ISFsiX) + 11F4(A)} 

+ 6A[M3], {F3(A) - F4(A)} + A[M%F4iX), (63) 

where [M"]c are the generalized susceptibilities for the initial condition satisfying 

= dzM{z)^ )■ 


III. EXTENSION TO MULTI-PARTICLE SYSTEMS 
A. Cumulants of net numbers 

Next, let us extend our model to a system composed of two species of Brownian particles. We then consider the 
cumulants of the difference of the particle numbers. This extension is needed for the description the cumulants of the 
net-baryon and net-electric charge nnmbers, which are defined by the difference of particle numbers. 

We denote the densities of two particle species as ni{x) and n 2 {x), and assume that their time evolutions are given 
by Eq. ED- We then consider the cumulants of the difference of the particle numbers in a range of width A, 

Q„et(A,t) =Qi(A,t)-Q2(A,t), (65) 

.A/2 

Qi{A,t) = I dxni{x,t). (66) 

J-A/2 

One easily finds that the distribution of Qnet approaches the Skellam one in the t —>■ oo limit, because in this limit 
the distributions of Qi and Q 2 are given by Poissonian and uncorrelated. 

For clarity now we recover the discretized notation for the particle density for the moment. We write the probability 
distribution function that a cell labeled by m have and rim^ particles at time t as Because rim^ 

and separately obey Eq. ED, if these two densities have no correlation in the initial condition, and are 
independent for any t. In this case, the probability can be factorized as 

P(n(i\n(2),t) = Pi(n(i),t)P2(n(2),t). (67) 

When is factorized in this way, the cumulant generating function K{ 6 i, 62 ,t) defined by 

0 ( 1 ) 0 ( 2 ) m 


is also decomposed as 


K{ei,e2,t) = Ki{ei,t) + K2{e2,t), 


(69) 


where K{6,t) is the cumulant generating function for the distribution fnnction of each particle. 

Using the generating function Eq. (1551) , cumulants of the difference of particle numbers hm = are given 

by 


c 


d^K 

dOmidOm^ ■ ■ ■ d9mi ’ 


(70) 
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where djdQ = dJdQx — i9/902- By combining Eqs. (l70t and (l6^ . one finds 


f. ...f. \ C_iy/^( 2 )^( 2 ) . . . „( 2 )\ 

"'■mi/c \ r; \'‘mi "-m2 "‘mil c- 


(71) 


In the continuum notation and after taking the integral over the range A, one obtains 


((g„et(A,t))"), =((gi(A,t))"), + (-i)"((Q 2 (A,t))”)„ 


(72) 


for initial conditions without correlation between two particle species. 

When the fluctuations of ni(x) and n 2 {x) are correlated in the initial condition, the probability distribution is no 
longer factorized as in Eq. (1671) . In this case, we denote the probability distribution as a superposition of the solutions 
of fixed initial conditions. 


P[ni{x),n 2 {x),t]= ^ F[Mi(z),M 2 (z)]PMi(z),M 2 (z)["i(a;),'" 2 (a;),t] 

{Mi(a:),M 2 (a;)} 

= X! F[Mi{z),M2{z)]PM^{z)[ni{x),t]PM^{z)[n2{x),t], (73) 

{Mi(x),M2(x)} 


where 7 mi( z), M 2 (z)[’2i(a;), 712(0;), t] are the solutions for the fixed initial condition with ni(x, 0) = Mi(x) and -"2(0;, 0) = 
M 2 (x). In the second line of Eq. (1731) we used the fact that ni (x) and n 2 (x) are uncorrelated with the fixed initial 
condition. It is also concluded from this uncorrelated nature that the cumulants of Pmi(z),M 2 (z) ["i(2;)i '" 2 ( 0 :), t] satisfy 
Eq. (17^ . Using Eq. dH]) and the relation for superposition of probabilities in appendix]^ the cumulants of gnet(A, t) 
are given by, 


(gnet(A,t))c = J dz{M^etiz))c,oHl{z), (74) 

((gnet(A,t))^)c = J dzidZ2{Mnet{zi)Mnet{z2))c,oHl{Zi)Hi{z2) 

+ J dz{Mtotiz))c,oH 2 iz), (75) 

((gnet(A,t))^)c = J dzidZ2dZ3{Mnet{zi)Mnet{z2)Mnet{z3))c,oHl{zi)Hi{z2)Hi{z3) 

Pi j dzidZ2{Mnet{zi)Mtot{z2))c,oHl{zi)H2{z2) 

+ J dz{Mnetiz))c,oH3iz), (76) 

r ^ 

((gnet(A,t))^)c = / dzidZ2dZ3dZ4{Mnetizi)Mnetiz2)Mnetiz3)Mnetiz4))c,oY\_Hi{zi) 

i—1 

+ 6 J dzidZ2dZ3{Mnet{zi)Mnet{z2)Mtotiz3))c,oHl{zi)Hi{z2)H2{z3) 

+ 3 J dzidz2{Mtot{zi)Mtot{z2))c,oH2{zi)H2{z2) 

+ ^J dzidZ2{Mnet{zi)Mr,et{z2))c,oHl{zi)H3{z2) 

+ J dz{Mtotiz))c,oH4iz), (77) 


with Mnet{z) = Mi{z) - M 2 {z) and Mtot{z) = Mi{z) + M 2 {z). 
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When the initial condition satisfies the locality condition Eq. (153, cumulants of (5net(A,t) are 

(Qnet(A,t))c=A[Mnet]c, (78) 

((Qnet(A,t))2), =A[Mtot]c{l - i"2(A)} + AJ,F2(X), (79) 

((Q„et(A, t)f), =A[M,et]c {1 - 3F2(X) + 2 F 3 (X)} + 3A[M„etMtot]c {F 2 iX) - FsiX)} 

+ A[M„yeE3(X), (80) 

((Qnet(A,<))4)e =A[Mtot]c{l " 7F2{X) + UF^iX) - 6F4(X)} 

+ 3A[(Mtot)"]c {i^ 2 (A) - 2 E 3 (X) + F 4 (X)} 

+ 4A[M„y, {F2(X) - 3F3(X) + 2F4(X)} 

+ {F3(X) - F4(X)} + A[M4,J,F4(X). (81) 

Here, [(Mnet)*(-^tot)'’]c are generalized susceptibilities defined similarly to Eq. (l58l) . 

The cumulants in equilibrium in this system is obtained by taking the large t limit in Eqs. ()78|) - (1811) . By taking 
X — > oo, Fn{X) for n> 2 vanish and the cumulants of Qnet become the Skellam ones with 

((Qnet)"*+')c = A[M„et]c, {{Qnetf^)c = A[Mtot]c, (82) 


with i being an integer. 


B. Generalization to multi-particle systems 

The above result would be suitable to describe the time evolution of net-baryon number cumulants in the hadronic 
medium in heavy ion collisions, because the net-baryon number, Xg®*, is given by the difference of baryon and 
anti-baryon numbers, Xb and Xg, respectively, as 

Xr' = XB-Xg. (83) 

When one considers the cumulants of net-electric charge, there are hadrons having ±2 charge in the unit of elemen¬ 
tary charge in addition to those having ±1. The net-electric charge is thus given by 

X^'^* = 2Xq,2 + Xq,i - Xq,_i - 2Xq,_2, (84) 

where Xq^„ is the number of hadrons with n electric charge. Owing to the existence of nonzero Nq^± 2 , then, even 
in the HRG model in which all Xq „ are given by the Poisson distribution and uncorrelated, the distribution of Xq®* 
deviates from the Skellam one. 

To describe this non-Skellam property, one may extend the DME to include four particle species. By constructing 
the net-electric charge as in Eq. (1841) , the distribution of the net-electric charge number in equilibrium has a deviation 
from the Skellam one. The effects of the charge 2 hadrons, however, are not large in the HRG model because of their 
small abundance [l^. We thus do not consider their effects in this paper and left the incorporation of these effects 
for future study. 


IV. RAPIDITY WINDOW DEPENDENCE OF CUMULANTS IN HEAVY ION COLLISIONS 

Next, using the result obtained in the previous section we consider the time evolution of the cumulants of conserved 
charges and their A?; dependence observed in relativistic heavy ion collisions. As for conserved charges, we consider 
the net-electric charge and net-baryon number, which are both observable in these experiments. Note that the 
measurement of the net-baryon number cumulants is possible BUS: although the detectors cannot observe neutral 
baryons. 

In QGD at sufficiently low temperature and low baryon chemical potential, the thermodynamic quantities including 
cumulants of conserved charges are well described by the HRG model. In this model, the distributions of net-baryon 
number and net-electric charge are given by the Skellam distribution to a good approximation. In dense and/or hot 
medium, the cumulants of conserved charges are expected to have deviations from the HRG values reflecting the 
change of the medium property including phase transitions to the deconfined medium. For example, fluctuations of 
net-baryon number and net-electric charge normalized by entropy density or net-baryon number are suppressed in 
the deconfined phase reflecting the decrease of the charges carried by quasi-particle excitations [iB- Near the QGD 
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critical point, the cumulants are expected to show characteristic suppressions and enhancements reflecting the critical 
phenomena [TJ, [l^, [11 ■ The goal of the experimental study of the cumulants of conserved charges is to And and 
confirm these signals [ 6 |. 

These signals, which would be developed in the deconflned medium or near the QCD critical point, however, are 
not directly observable. Before they are observed by the detectors, the hot medium undergoes a time evolution 
in the hadronic stage. The fluctuations are modified in the hadronic medium, presumably toward the one of the 
equilibrated hadronic medium. When one interprets the experimental results on fluctuations, one must keep this 
effect in mind; direct comparisons of the experimental results on cumulants with a thermal value at some early stage 
in heavy ion collisions would lead to a wrong conclusion. The purpose of this section is to describe the diffusive 
process of fluctuations in the hadronic stage and to estimate their effects on experimental results using the solution 
of the DME obtained in previous sections. 

In relativistic heavy ion collisions with sufficiently large collision energy per nucleon, .^snNi the hot medium 
created by the collisions has an approximate boost invariance around mid-rapidity. Useful coordi nates to describe 
such a system are the coordinate-space rapidity y = (1/2) log(t -I- x){t — x) and proper time r = — x'^, where x 

denotes the longitudinal direction. In this coordinate, the diffusion equation in the Cartesian coordinate, Eq. @, is 
rewritten as 


d 

—n(77,r) = Dy{T)—n{y,T), 


(85) 


where n{y,T) is the density per unit rapidity, and the diffusion constant in this coordinate, Dyir)^ is related to the 
Cartesian one D as 


Dy{T)^DT-\ ( 86 ) 

Because of Eq. (1551) . the diffusion constant D has to be replaced by Dy when one applies the result in the previous 
sections to the description of the diffusive process along the rapidity direction. 

In experiments, the fluctuations of particle numbers in a pseudo-rapidity window A 77 are observed after integrating 
out the transverse direction ii- The pseudo-rapidity window can be varied within the coverage of the detector. 
In the following, we regard the pseudo-rapidity window as the one of the coordinate space rapidity rj [s^. The 
A 77 dependence of the cumulants can then be compared with the A dependence obtained in previous sections. The 
cumulants obtained in the previous section, Eqs. (|78| - dn]), depends on the initial condition and the diffusion distance 
d{t). By fitting the experimental data by these results, one can constrain the initial condition and the value of d(t), 
namely information on the thermodynamics and transport property of the hot medium. 

In heavy ion collisions, the conserved charges in the total system is fixed and does not fluctuate. Because the hot 
medium created in heavy ion collisions is a finite volume system, the measurement of fluctuation in a sub-volume is 
affected by this effect. Eor sufficiently large .^/snn and moderate range of A? 7 , however, such an effect is suppressed 
because of the finite duration of the diffusive process in the hadronic medium [^ . 


A. Initial condition 

In this study, we set the medium at chemical freezeout time, r = tq, as the initial condition, and apply the 
DME to describe the diffusive process in the hadronic medium until kinetic freezeout. Eor sufficiently large 
hadronization and chemical freezeout are expected to take place at almost the same time. Because the local density 
profile of conserved charges in each event is almost frozen within the short duration because of the charge conservation, 
if the fluctuations are well equilibrated in the early stage, the fluctuations of conserved charges at r = tq will strongly 
reflect the thermal property in the deconflned medium. The fluctuations of non-conserving quantities, on the other 
hand, are not constrained by the conservation law, and thus would be sensitive to the hadronization mechanism rather 
than the primordial thermodynamics. 

In this study, we assume that the locality condition Eq. dSZl) is satisfied at t = tq. Note that this assumption is 
satisfied if the medium is fully equilibrated at this time. With this assumption, the initial condition is specified by the 
generalized susceptibilities. To describe the cumulants up to fourth order, we need six susceptibilities in the initial 
condition 

[MLtU [MLU (87) 

The first three susceptibilities are theose of conserved charges. By normalizing them by their equilibrated values given 
in Eq. (15^ . we introduce the following three parameters 

D2 = [M^^^]c/[Mtot]ct D3 = [M)(g^]c/[Mnet]c, D4 = [M^^^]c/[Mtot]c, 


( 88 ) 
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Ar|/d(x) ATi/d(T) 

FIG. 3. Ari/d{T) dependences of the functions defined in Eqs. (1961) - (19911 . 


representing the magnitude of the cumulants of conserved charges at hadronization. Other three parameters are not 
the cumulants of conserved charges. We parametrize these quantities by 

a = [Mnet.^tot]c/[.Wnet]c, b = [M^^^Mtotlc/ [Mtot]c, C = [(Mtot)^]c/[.Mtot]c, (89) 

where we again normalize them by the equilibrated values. Among these parameters, Dn are the quantities which 
should be compared with theoretical studies, such as lattice simulations, on the cumulants of conserved charges. Note 
that D 2 and c are positive definite because they are second-order fluctuations, while other parameters can take both 
positive and negative values. In addition to Eqs. (IMll and (IMll . the A 77 dependence of the cumulants depends on the 
the diffusion length d(r). 

Because the parameters in Eq. (IMl) are not directly constrained by the conservation law, their values would be 
sensitive to the hadronization mechanism rather than the primordial thermodynamics. In particular, the parameter c 
does not include the net-particle number and is completely free from the conservation law. Experimental constraint 
on this parameter would provide us with novel information on the hadronization mechanism |19| . 

Eor net-baryon number cumulants and for sufficiently low energy collisions at which the creation of the anti-baryons 
are well suppressed, the system can be regarded as the single particle one without anti-baryons. Because M^et and 
Mtot become identical in this case, one obtains a = c = D 2 and b = D 3 . The number of parameters are reduced in 
this case. 


B. Normalized cumulants 


In the following, we plot the Ary dependences of cumulants by normalizing by their equilibrated values. 


Rn{X) = 


((Qnet(A?y,T))")e 
limr-i.oo((Qnet(A?y, r))”)c 


(90) 


We refer to Eq. (1^01) as the normalized cumulants in what follows. Note that the equilibrated values are the Skellam 
one given in Eq. (1821) in our model. 

When one analyzes the normalized cumulants in experiments, one may use the equilibrated values in the HRG 
model for normalization. An alternative way to deduce the equilibrated value is to take the Ary —>■ 0 limit of the 
experimental results on the cumulant. 

Usually, the experimental results on the higher order cumulants are discussed in terms of their ratios such as R 4 /R 2 
and R 3 /R 1 0. When one considers the Ary dependences, however, the normalized cumulants would be better than 
the ratios, since one can investigate the individual cumulants directly using i?„. 

Using the parameters in Eqs. (1551) and (1591) . the solution obtained in the previous section Eqs. (EHl) - dH]) in terms 
of the normalized cumulants are given by 


R2{X) = G2fi{X) + £>2^2,i(A), 

RsiX) = G3 ^o{X) + D^Gs^ilX) + aG3_2(A), 

RilX) = G4 _q{X) + D^G/i^ilX) -|- D2Gi^2{X) -\- bG^^^^X) + cG^^^lX), 


(91) 

(92) 

(93) 
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FIG. 4. Ari/d{T) 


1/X dependence of the normalized second order cumulant R 2 for D 2 — 0 and 0.5. 


where we have introduced the following functions: 


G2 ,o{X) = 1-F2{X), G2,i{X) = F2 {X), (94) 

G3,o(^) = 1 - 3F2(X) + 2F3iX), GsAX) = F^iX), (95) 

GsMX) = 3(F2(X) - F3(X)), (96) 

Gi,o{X) = l-7F2iX) + l2F3iX)-6F4X), G^^X) = Fi{X), (97) 

Gi^X) = 4(F2(X) - 3F3(X) + 2F4X)), (98) 

G4,3(^) =6(F3(X)-F4(X)), G4,4(X) = 3(F2(X)-2F3(X)+F4(X)). (99) 


In Fig. [U we show 1/X = Ari/d{T) dependence of the functions in Eqs. (|M)l - (IM|l . The figure shows that these 
functions behave differently as functions of Arj/d{T). The structures of these functions are responsible for characteristic 
behaviors of Ar] dependences of the cumulants discussed in the next subsections. 


C. Second order cumulant and diffusion distance 

Now, let us examine Ar]/d(T) = 1 /X dependence of the normalized cumulants Eq. (lOTl . We first consider the the 
second order one i? 2 - As in Eq. m, i ?2 depends on the initial condition only through Z? 2 . In Fig. 0] we show the 
Ar]/d{T) dependence of R2 for D2 = 0 and 0.5. 

Since we plot the Aij dependences of the cumulants as functions of Arj/d^r) throughout this section, it is instructive 
to give a rough estimate on the value of ^(t) at this point by comparing Fig. 3] with the existing experimental results. 
The second order cumulant of net-electric charge has been observed by ALICE collaboration @ . In Ref. @ , the result 
is plotted using the quantity called the D-measure @, which is related to R 2 as 

D = AR2, (100) 

provided that Qnet is the net-electric charge. IW comparing the result in Fig. H] with Fig. 3 in Ref. one can 
constrain the values of D 2 and (i(r) [2l|. In Ref. the value of £>2 is estimated as D 2 = 0.5. Using this value of D 2 
and Fig. 3 in Ref. @ one can estimate [2l| 

d(r) = 0.3 - 0.5. (101) 

This result shows that the maximum rapidity window of the ALICE detector, A?7niax = 1-6 corresponds to 

A??„,ax/d(r) = 3.2 ^ 5.3. (102) 

The ALICE detector thus can analyze Ar]/d{T) dependences of cumulants for 0 < Ar] < Arjn^ax- Note that this 
estimate on (i(r) strongly depends on the value of £> 2 ; for smaller value of D 2 , d^r) becomes much larger. 

When one applies the net-baryon number to Qnet, the value of of dir) should be much smaller because the hadrons 
having baryon number, namely baryons, are considerably heavier than those having electric charge, which are dom¬ 
inated by pions in heavy ion collisions. This suggests that wider range of Arj/dir) can be analyzed for net-baryon 






17 
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FIG. 5. Normalized third-order cumulants R 3 as a function of 1/X 


Arj/d{T) for D 3 = —0.5, 0, 0.5 and several values of a. 
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FIG. 6. Normalized fourth-order cumulants R 4 as a function of Ar)/d(T) for various initial parameters D 4 , D 2 , b and c. The 
left (right) panel shows the result for D 2 = 0 {D 2 = 0.5). 


number cumulants with a fixed A 77 coverage of an experimental detector. The use of the net-baryon number cumu¬ 
lants is advantageous compared with net-electric charge in this sense. The simultaneous measurements of both the 
net-baryon number and net-electric charge cumulants in a same experiment would provide us interesting results, such 
as the difference of the diffusion constants of net-baryon number and net-electric charge. 


D. At; dependence of higher order cumulants 

1. Initial condition with small fluctuations 

Next, we focus on the Arj dependences of third and fourth order cumulants. We first consider the initial condition 
with small fluctuations, since the fluctuations of net-baryon and net-electric charge numbers are expected to be 
suppressed in the deconfined medium The higher order cumulants can also be suppressed near the QCD 

critical point reflecting the critical phenomena 

The normalized cumulant R 3 depends on ZI 3 and a, while R 4 depends on D 4 , D 2 , b and c. In Fig. O we show 
dependence of R 3 on Ari/d{T) = 1/X with D 3 = —0.5, 0 and 0.5 for several values of a. In Fig. [HI A-q/dij) dependence 
of i ?4 is shown for various values of the initial parameters. The figures show that the dependences of i?„ on Arj/dir) 
are sensitive to the initial parameters in the range Ar]/d{T) < 7. For example, one sees from Fig.[S]that the behavior 
of i ?3 is sensitive to a for Ar]/d{T) < 1 while the effect of D 3 becomes dominant as Ar]/d{T) becomes larger. This 
behavior indicates that the experimental measurement of R 3 in the range Arj/d{T) ^ 5 makes it possible to constrain 
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FIG. 7. Dependences of the normalized non-Gaussian cumulants R 3 (left) and R 4 (right) on /S.r]/d(r) with D 2 = D 3 = D 4 — 1. 


two parameters ZI 3 and a. It is also notable that Arj dependence of R 3 can become non-monotonic for several choices 
of parameters. Experimental observation of such a non-monotonic behavior should be an interesting and unique signal 
to constrain these parameters and underlying physics. Figure [5] also shows that the behavior of R 4 is sensitive to the 
initial parameters. For example, R 4 can have two extrema for a choice of initial parameter in the range Arj/dir) < 5. 
Although there are 16 lines for R 4 in each panel of Fig. IH these lines would be experimentally distinguished if the 
experimental result on R 4 is obtained within 10% precision. This precision will be achieved in the BES-II program 
at RHIC [I|. Such experimental analysis can constrain the initial parameters D 4 , b and c as well as ^(t). We also 
note that, although R 4 depends on D 2 ^ the value of D 2 can be constrained by the measurement of i? 2 - Similarly, 
because all normalized cumulants i? 2 , R 3 and R 4 depend on a common d{T), this parameter can be constrained by 
the simultaneous use of these three cumulants. In this way, the complementary use of these cumulants will lead us 
to a deeper understanding of the experimental results. The analyses of both the net-baryon and net-electric charge 
fluctuations in a same set of collision events will also give us further interesting information on the difference of the 
diffusion constants for these charges and so forth. 

For the analysis of the net-baryon number cumulants with small because the abundance of anti-baryons are 

suppressed the parameters are constrained as a ~ D 2 , c ~ D 2 and 5 ~ Ua. Using these conditions, the analysis of the 
Ar] dependence of normalized cumulants i? 2 , R 3 and R 4 would further constrain the physics of fluctuations. 

It is worth emphasizing that the results in Figs. [5] and | 6 ] show that the value of a cumulant with a fixed Arj can 
deviate from both the initial and equilibrated values significantly. In particular, the sign of non-Gaussian cumulants 
can become negative even when the initial and equilibrated values are positive. These results show that the cumulants 
observed in experiments with a fixed A 77 should not be compared with the theoretical analysis of the cumulants 
obtained in statistical mechanics assuming an equilibration. Total description of the Arj dependence of the cumulants 
are indispensable for proper understanding of the experimental results. 


2. Initial condition with large cumulant 

Next, we show the Arj dependence of the cumulants for initial conditions where higher order cumulants of conserved 
charges are consistent with or larger than the equilibrated values. 

We first consider the case D 2 = = D 4 = 1, i.e. the cumulants of the conserved charge in the initial condition is 

given by the Skellam distribution. In Fig. [3 we show the At; dependences of third and fourth order cumulants with 
this initial condition. In the figure, we vary the parameters a, b and c, which are not constrained by the conservation 
law. The figure shows that the normalized cumulants for a given Ar] can have a significant deviation from unity 
depending on the initial parameters although initial and equilibrated values of the cumulants are both unity. Only 
when a = 6 = c = 1, all normalized cumulants become unity irrespective of the value of A 77 . This result again tells us 
that the measurement of the cumulants with a fixed Arj should not be compared with those in equilibrated medium. 

Finally, we consider the At; dependence of the normalized cumulants R 3 and R 4 for relatively large D 3 and D 4 . 
Such a case would be realized when the medium near the critical point is formed in the time evolution of the hot 
medium and the critical enhancement of fluctuations are well developed m- In Fig. |S]we show the At; dependences 
of i ?3 and i ?4 with D 3 = 4 and D 4 = 4, respectively, for several values of initial parameters. The figure tells us the 
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ATi/d(T) Ar|/d(x) 

FIG. 8. Dependences of the normalized non-Gaussian cumulants R3 (left) and R4, (right) on Arj/d{T) with D3 = 4 and D4 = 4 
with various initial parameters. 


same conclusions as before: (1) The values of R 3 and R 4 with a fixed Arj can have significant deviation from the 
initial and equilibrated values, and (2) initial parameters can be constrained experimentally from the Arj dependences 
of the normalized cumulants. 


V. DISCUSSIONS AND SUMMARY 

In the previous section we have discussed the rapidity window Arj dependence of the cumulants. These results can 
be directly compared with the experimental results of the cumulants of net-baryon number and net-electric charge. 
By this comparison, especially by the simultaneous use of the three normalized cumulants R 2 , R 3 and i? 4 , parameters 
in the model and underlying physics on the thermal and transport properties of the hot medium can be constrained 
experimentally. The discussion in the previous section also tells us the importance of the total understanding of 
the A 77 dependences of the experimental results. Because the value of the cumulants with a fixed Arj can deviate 
significantly from the initial and final values, they should not be directly compared with theoretical studies which are 
obtained in statistical mechanics assuming equilibration, such as lattice simulations. The quantities which should be 
compared with the cumulants in statistical mechanics are D 2 , D 3 and D 4 , which can be deduced by the analysis of 
the A 77 dependence. 

We note that the Ary dependence of the Gaussian fluctuation encodes the same physical information as those 
accessible with the balance function [42l - l43| of the corresponding particle species, because the former is obtained 
by the integral of the latter, and vice versa [U . On the other hand, the Ary dependences of non-Gaussian 

cumulants are related to multi-particle correlation functions, and provide us with a novel information on the diffusive 
process of conserved charges. 

All the results in this study is obtained using the DME. In this model it is assumed that the Brownian particles in 
the model do not interact with each other. This assumption would be qualitatively justified when Qnet is identified to 
be the net-baryon number because of the following two reasons. First, baryons in the hadronic medium mainly interact 
with pions, and the baryon-baryon interactions scarcely take place [l7|. Second, the success of the thermal model for 
the particle abundances indicates that the pair annihilation of baryons is well suppressed after chemical freezeout [II- 
The baryons thus can be regarded as the Brownian particles floating in the pionic medium without interacting with 
one another, and the correlation between two baryons in the hadronic medium will be well suppressed. For net-electric 
charges, on the other hand, the interactions, especially the pair-annihilation, occur frequently between hadrons having 
electric charge, especially pions. The present model thus would not be suitable for the quantitative description of 
the diffusive process of the net-electric charge fluctuations. It is desirable to extend the model to incorporate the 
interaction between particles for more quantitative description of the net-electric charge fluctuations. In the present 
study it is also assumed that the fluctuations satisfy the locality condition at chemical freezeout time, which, however, 
might be violated in realistic collision events. To incorporate such effects one has to modify the initial condition. We 
left these extensions to future study and do not address further in this study. 

In this paper, we considered the time evolution and rapidity window, A?y, dependences of the cumulants of conserved 
charges observed in relativistic heavy ion collisions. In order to describe the non-equilibrium time evolution of non- 
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Gaussian cumulants we employed the diffusion master equation Eq. O- The analytic formula for the time evolution 
of cumulants with arbitrary initial conditions are obtained. Using these results, we discussed the Arj dependence of 
cumulants of conserved charges observed in heavy ion collisions. Various suggestions have been made to utilize the 
Ar] dependence of the non-Gaussian cumulants for the understanding of the thermal and transport properties of the 
hot medium. 

The author thanks Masayuki Asakawa, Hirosato Ono and Miki Sakaida for valuable discussions. He also thanks 
Shinichi Esumi, Volker Koch, Hiroshi Masui and Krzysztof Redlich for useful conversations and encouragements. This 
work was supported by JSPS KAKENHI (Grant Number 25800148). 


Appendix A: Superposition of cumulants 


In this appendix, we consider cumulants of a probability distribution function which is given by a superposition of 
probability distribution functions. Let us consider a probability distribution function P{x) for an integer stochastic 
variable x, and assume that P{x) consists of the superposition of sub-probabilities as 

P{x)=Y,F{N)PNix), (Al) 

N 

where Pn{x) are sub-probabilities labeled by integer N. Each sub-probability is summed with a weight F{N) satisfying 
F{N) = 1 , which is also regarded as a probability. The purpose of this appendix is to represent the cumulants 
of P{x) using those of Pn{x) and F{N). Although we write down the results explicitly up to fourth order in this 
manuscript, the result can be extended to higher orders straightforwardly. A similar result given in this appendix 
is presented and used in Ref. [13 to relate the net-proton and net-baryon number cumulants in the final state in 
relativistic heavy ion collisions. 

We start from the cumulant generating function of Eq. m, 


K{9) = log ^ e^^P{x) = log XI H (^2) 

X N X 

= log ■£nN)Y: f.KN{e) (AS) 

N X 


where Kn{ 9) — \og^^ e^^PN{x) is the cumulant generating function for Pn{x). Using the cumulant expansion, 
Eq. (jAS|) is written as 


K{e)=J2^Y.i^N{e)]T (A4) 

m F 

= J2Kn{9) + i Y.^SKM{9)y + i 

F F F 

+ ^ X[AlAf(^)]e 4 I (A5) 

F 


where J2 f 1® 1"^® shorthand notation of J2 n 1® 1^^® mth-order cumulant of ATjv for the sum 

over F, whose explicit forms up to fourth order are given in the far right hand side with 

F F V F / 

= Y.(SKN{9) f - 3 ■ (A7) 

F F \ F / 

Cumulants of P{x) are given by derivatives of K{9) as 

{x-)e = ^K{0) = K(-\ (A8) 

All cumulants can be obtained with Eqs. (IA8|) and (IA5I) . In order to calculate the cumulants explicitly, we first note 
that the normalization condition PnIx) = 1 yields Kn{9) = 0. From this property, it is immediately concluded 


that all Kn{9) in a term in the far right hand side of Eq. (lASp must receive at least one differentiation so that the 
term gives nonzero contribution to Eq. (IA 8 I) . This means that the mth order term in Eq. (IA5I) can affect Eq. (IA 8 I) 
only if m < n. Keeping this rule in mind, derivatives of Eq. (IA5I1 with 0 = 0 is given by 


(A9) 

F 

=Y,K^n^+Y.^5K^n)\ (AlO) 

F F 

K(3) = Y,Kf (All) 

F F F 

“ E «J;’"4”+3 +6 

F F F F 

+ (A12) 

F 

with = d^KN{0)/d9^ being the cumulants of sub-probabilities Pn{x). Equations (IA9I) - (IA12I1 relate the 

cumulants with 

The above relations are further simplified when the cumulants of Pn{x) are at most linear with respect to N, i.e. 

=Are(„)+C(n), (A13) 

where ^(„) and C(n) are constants which do not depend on N. Substituting Eq. (IA13I1 in Eqs. (IA9I1 - (IA12I) one obtains 

=C(i) +C(i)(-A)f, (A14) 

=C(2) + ii2){N)F + (A15) 

=C(3) + C(3)(-A)f + H(i)£,(2){^N‘^) F + ^(i)^{5N^)f, (A16) 

AT^A _|_ ^(4)(7V)^ -I- (4^(i)^(3) -I- 3^(2)^)(5A^^)_F -I- £,(2){5N^)F 

+ (A17) 

where {0{N))f = Sat 0(-A)A(A^) denotes the average over F{N); these averages in the above formulas represent 
cumulants of the probability F{N). 

These results are easily generalized to the case where the sub-probabilities are labeled by multi-dimensional integers, 

P{x) = J2FiN)PN{x), (A18) 

N 

with N = (A^i, A^ 2 , * * * , ^m) representing a set of m integers. When the cumulants of Pn{x) are linear in i.e., 

=iV •€(„)+ C(n), (A19) 

a similar manipulation as before yields, 

=({!)+{iwN)F, (A20) 

=C(2) + (1(2) • N)f + ((Id) ■ SN)^)f, (A21) 

F^^ =C(3) + (1(3) • N)f + 3((|(1) • <5iV)(^(2) • SN))f + ((Id) • SN)^)f, (A22) 

kW =C(4) + (Id) • N)f + 4((|(i) • <5iV)(|(3) • SN))f + 3((|(2) • 6N)^)f 

+ 6((|(i) • <5iV)2(|(2) • 6N)}f + ((Id) • SN)^),,f. (A23) 

Appendix B: Chemical Reaction 


In this appendix, we consider a simple chemical reaction between two species of particles X and A [3^ . 
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It is assumed that a particle X (A) becomes A (X) with a probability ki (^ 2 ) per unit time. We denote the number 
of each particle in the system as x and a. For simplicity, it is further assumed that a is sufficiently large compared 
with X so that the number a can be regarded fixed. The time evolution of the probability distribution P{x, t) is then 
described by the master equation [s^ 

dtP{x, t) = k 2 aP{x — 1, t) + ki{x + l)P{x + 1, t) — {kix + k 2 a)P{x, t). (B2) 

In the following, we present the solution of this master equation for arbitrary initial conditions. 


1. Fixed initial condition 

Let us first solve Eq. (IB2() with a fixed initial condition 

Pix,0)=S,,N, (B3) 

i.e. the number x is fixed to N without fluctuations at t = 0. 

In order to solve Eq. (IB2|) . we use factorial generating function, 

Gf(s,f)=^s"P(x,t). (B4) 

X 

Substituting Eq. (IB4I) in Eq. (IB2I) . one obtains 

dtGf{s, t) = k 2 a{s — l)Gf (s, t) — fci(s — I)9sGf (s, t). (B5) 

Equation (IB5I) can be solved by using the method of characteristics, which yields 

Gf (s, t) = F((s - (B 6 ) 

with an arbitrary function F and the average number of x in equilibrium, Neq = k 2 a/ki. 

The function F is determined by specifying the initial condition. The generating function corresponding to the 
initial condition Eq. (IB3I) is Gt(s,0) = . This gives 

E(r) = (r + I)^e’'"‘='>. (B7) 

The solution of Eq. (IB5I) with the initial condition Eq. (IB3I) thus is 

Gf(s,t) = (l + (s- (B 8 ) 

The moment and cumulant generating functions, G{6,t) and K{6,t), respectively, are given by 

G{B,t) = Gi{e\t), (B9) 

K{e,t) =\ogG{e,t). (BIO) 

Derivatives of K{9,t) define cumulants; {x{t)^)c = K{0, t)/. Using 

^0(9, () = (), (B12) 

and so forth, one obtains 

{x{t)), = A^eq(l - (B13) 

{x{tf), = fVeq(l - (B14) 

{x{tf)^ = A^eq(l - 6“'=^*) + A^(e-'=1‘ - 3e-2'=i* + 2e-3'=i*), (B15) 

{x{tf)c = iVeq(l - 6“'=^*) + N{e-^^^ - 7e-2fci* + 126-3'=^* - (BI6) 
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2. General initial condition 

In order to obtain the t dependence of the cumulants for general initial conditions, we make use of the result in 
Appendix El Since Eq. (IB2I) is a linear differential equation, the solution of P{x,t) for an initial condition P{x,0) = 
F{x) is given by 


Pix,t)=Y,FiN)PNix,t), (B17) 

N 

where Pn{x, t) is the solution for the fixed initial condition, Eq. (IB3I) . We further remark that the cumulants Eqs. (IB13I) 
- (jB16|) are at most linear in N and fulfill the form in Eq. (jA13p with 

C(„) = fVeq(l - (B18) 

and etc. The cumulants for arbitrary initial conditions thus are obtained by substituting Eqs. (IB18|) in Eqs. (IA14I) - 
(I A 171) . The results up to fourth order are 

{x{t)), =7Veq(l - (B19) 

=7Veq(l - + (IV)oe-'=^‘ + {{5N^)^ - (B20) 

{x{tf), =7Veq(l - 

+ ((<5iV3)o - 3(<5iV2)o + 2(7V)o)e-3'=i‘, (B21) 

(x(t)4), =7Veq(l - + (IV)oe-'=i‘ + - {N)o)e-^^^* 

+ 6{{6N^)o - 3{SN^)o + 2(7V)o)e-3'=i‘ 

+ {{SN^)c,o - 6{SN\ + n{SN^)o - 6(iV)o)e-4'=i‘, (B22) 

where {0{N))q denotes the average in the initial condition. 

From Eqs. (IB19I) - (IB22I) . one finds that in the large t limit, kit —?> oo, all cumulants converge a same value, 

(cr(t)”)e = iVeq, (B23) 


which means that the distribution approaches the Poissonian in the large t limit irrespective of the initial condition. 


3. Some examples of time evolution 

It is instructive to see the time evolution of {x(t)'^)c with some specific initial conditions. 


a. Poisson distribution 

First, let us consider the initial condition that P{x,t) is the Poisson distribution at t = 0. In this case, the initial 
condition satisfies 


{N)o = {SN^)o = {SN\ = (5iV4),,o. (B24) 

Substituting Eq. (IB24I) in Eq. (IB19I) - (IB22I) , all terms in the right-hand side vanish except for the first two terms and 
one finds that the all cumulants have the same t dependence 

{x{tr), = IVeq(I - + (IV)oe-'=^‘. (B25) 

This result shows that the distribution P{x,t) starting from the Poissonian stays Poissonian for all t [s^, while the 
average of the distribution shifts from {N)o at t = 0 to iVeq for t —^ oo. 


b. {N)o = iVeq 


Next, we consider the initial condition 


P{X,0) = 4,Weq, 


(B26) 
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FIG. 9. Dependences of the cumulants {x{t)'^)c on t with fixed initial condition Eq. (IB26II . 


i.e. the particle number is fixed to the equilibrated value Ngq without fluctuation. In this case, one obtains 


{x{t))c =Neq, 

(B27) 

{x{tf)c =iVeq(l - 

(B28) 

{x{tf)c =iVeq{l - + 2e-3'=i‘), 

(B29) 

(x(t)4)c =fVeq(l - + 12e"3'=^‘ - 66"^'=^*). 

(B30) 


In Fig. [3 we show the t dependence of second to fourth order cumulants for the initial condition Eq. (IB26I) up to 
fourth order. The figure shows that all cumulants approach N^q as t increases. One also finds that the growth of 
higher order cumulants is slower. This behavior is understood as follows. First, the higher order cumulants are more 
sensitive to the tails, i.e. far side from the average value, of the distribution. Second, propagation of the probability 
distribution to far side from the Ngq should take longer time. Thus, the approach of higher order cumulants to the 
equilibrated value is much slower than the lower order ones. 


Appendix C: Ix{C) and Fx 


In this appendix, we summarize property of Ix{C) defined in Eq. (HD) 


. 1/2 

IxiC) = / 

J-Xjl 

and Fn{X) defined in Eq. (ICT|) . 

By first integrating out p in Eq. m one obtains 

.1/2 


27T 


IxiC)= I 


1 - 1/2 
= - (erf 


'V^x 

iC+l/2 X-l/2 


V2X 


— erf- 


V2X 


Eor X = 0, IxiC) satisfies 


/■1/2 r fih , , /•1/2 

^o(C) = / dC = / dCSiC - C) = 0(1 - C") 

J-l/2 J 27r J-l/2 


where 0(a;) is the step function. From this result, one obtains 

[7o(c)r = 0(i-c"), 


(Cl) 


(C2) 


(C3) 


(C4) 
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FIG. 10. Functions F„{X) with n = 2, 3 and 4. 


for arbitrary n. 

In Sec. m we define 


From Eq. (El, one easily obtains 


/ OO 

dz[Ixiz/A)r. 

-OO 


FniO) = 1 


for any n. The integral for n = 2 can be performed analytically and one obtains 

For n = 1, Fi(X) is a constant 


Fi{X) = l. 


(C5) 


(C6) 


(C7) 


(C8) 


For n> 2, Fn{X) are monotonically decreasing functions of X with F„(0) = 1 and limx-i.oo Fn{X) = 0. In Fig. (TUI 
we show the X dependences of F„(X) for n = 2, 3 and 4. 
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